

/*
 *  main0.cpp
 *  
 *
 *  Created by Camilo Ruiz Méndez on 29/11/08.
 *  Copyright 2008 __MyCompanyName__. All rights reserved.
 *
 */






#include <iostream.h>
#include <complex.h>
//#include <fftw.h>
#include <alloc.h>

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define complex complex<double>
#include <fstream.h>

#include "../include/constants.h"
#include "../include/laserfield.h"
#include "../include/gridobject.h"
#include "../include/carpetas.h"
#include "../include/routines.h"
using namespace std;
void interpp(double *a);


int main()
{
	
	/*************************************************/
	//Create folders and files.
	
	FILE *fp0,*fp1,*fp2,*fp3,*fp4,*fp5,*fp6,*fp7,*fp8,*fp9,*fp10; //Fichero de prueba
	int dir_limit=50; //Numero de maximo de directorios
	
	//Creating a new data directory
	
	if(CrearCarpeta(dir_limit)==0)
	{	
		printf("Has superado el limite de directorios que tu pusiste, cazurro!!");
		exit(1);
	}
	
	fp0=fopen("datapic.txt", "wt+");
	if(fp0==NULL)
	{
		printf("Error opening the file.\n");
		exit(2);
	}
	fp1=fopen("datapic1.txt", "wt+");
	fp2=fopen("datapic2.txt", "wt+");
	//fp2=fopen("datapic2.txt", "wt+");
	fp3=fopen("datapic3.txt", "wt+");
	fp4=fopen("datapic4.txt", "wt+");
	fp5=fopen("datapic5.txt", "wt+");
	fp6=fopen("datapic6.txt", "wt+"); //Fichero creado por carlos para los primeros tests de armonicos
	fp7=fopen("datapic7.txt", "wt+"); //Fichero creado por carlos para los primeros tests de armonicos
	fp8=fopen("datapic8.txt", "wt+"); //Fichero creado por carlos para los primeros tests de armonicos
	fp9=fopen("datapic9.txt", "wt+"); //Fichero creado por carlos para los primeros tests de armonicos
	fp10=fopen("datapic10.txt", "wt+"); //Fichero creado por carlos para los primeros tests de armonicos
	
	/*************************************************/
	/* Define variables
	/*************************************************/
	
	
	int iaux,npion,npim;
	
	int np=50;//400;//3; //4000;              ///< Particle number
	int ndg=1024; //200; //1024*8;           ///< Number of points in the spatial grid
	int ndt=1024;//*4; //5000;//1024*32;      ///< Number of points in the temporal grid
	int nplim=2;

	
	double ran1;
	int iff,iseed;
	
	//Physicis variables
	double qm,ae;
	
	//Material parameters:
	double n0,nc,wplasma;
	double xgleft,xgright,lg;
	int ileft,iright;
	
	//Field parameters
	double E0,w0,k0,wp,kp,lambda,E0x,E0z,B0x,B0y,B0z;
	int nw;
	
	//Particle parameters
	double qp[np],mp[np];
	int na;
	
	//Ionization parameters
	double coef, exponente,cargaminima;
	int npminimo;
	
	//Constants
	//double pi=4.*atan(1.0);
	//complex iu=complex(0.,1.);
	//double unsexto=1./6.;
	//double untercio=1./3.;
	//double c=137.037;
	//double q=1.;
	//double mass=1.;
	
	int iflagion=0;
	int labelx=1, labely=2, labelz=3; //labels used as arguments for the function CalculaJ
	
	//Field parameters
	w0=0.057;
	k0=w0/lightC_au;
	lambda=2.*pi/w0;
	nw=20;
	
	wp=w0/2./nw; //Envelope frequency;
	kp=wp/lightC_au;
	E0=50.5;
	
	laser laser1;
	laser1.initialize(w0,nw,E0);
	
	//Calculation parameters.
	lg=(nw+1)*lambda/4; //WATCH Fix here to make deltag=c*deltat
	double tmax=(nw+1)*2*pi/w0;
	double deltag=lg/ndg;
	double deltat=tmax/ndt;
	
	grid space;
	grid time;
	space.set_grid(ndg,deltag);
	time.set_grid(ndt,deltat);
	
	
	//Material parameters
	xgleft=lg/2.-lambda/2.;
	xgright=xgleft+lambda/1.;
	nc=w0*w0*mass/4./pi/q/q;  //Volumetric density
	n0=7*nc;					//Volumetric density
	
	//Limits of the material sheet
	ileft=int(xgleft/deltag)+1;
	iright=int(xgright/deltag)+1;
	
	int if1=ileft;
	int if5=iright;
	int if3=int(0.5*lg/deltag)+1;
	int if2=(if1+if3)/2;
	int if4=(if3+if5)/2;
	na=3;
	qm=-q/mass;
	ae=0.5*qm*deltat;
	
	
	//complex IU;
	
	//Variables for output
	int nscleft, nscright;
	//int if1,if2,if3,if4,if5;
	int nsc;
	int idj[10000];
	
	//Arrays over the spatial grid.
	grid_object oRo(space), oRop(space),oJx(space),oJy(space),oJz(space), oRoold(space), oRopold(space);
    grid_object J[3];
    
	double Ro[ndg], Rop[ndg],Jx[ndg],Jy[ndg],Jz[ndg], Roold[ndg],Ropold[ndg];
	double p0[ndg],deltap0[ndg],wrate[ndg];
	double IntRo[ndg], IPhiM[ndg], IPhiP[ndg];
	double IntJx[ndg], IntJy[ndg], IntJz[ndg];
	double AuxIntRo,AuxIntJx, AuxIntJy,AuxIntJz;
	
	double IAxM[ndg], IAxP[ndg],IAyM[ndg],IAyP[ndg],IAzM[ndg],IAzP[ndg];
	double Ax[ndg],Ay[ndg],Az[ndg],phi[ndg];
	double DAx[ndg],DAy[ndg],DAz[ndg];
	double Axold[ndg],Ayold[ndg],Azold[ndg];
	double Exold[ndg],Eyold[ndg],Ezold[ndg],Byold[ndg],Bzold[ndg];
	
	double Ex[ndg],Ey[ndg],Ez[ndg],By[ndg],Bz[ndg];
	
	double t,xg;//,deltag,deltat,gamma;
	double xp[np],vpx[np],vpy[np],vpz[np];
	double xpnew[np],vpxnew[np],vpynew[np],vpznew[np];
	double vvx,vvy,gvxs,gvys,vvz,gvzs,dcg,fb;
	double aexpt,aeypt,aezpt,abxpt,abypt,abzpt;
	double deltavx[np],deltavy[np],deltavz[np],deltax[np];

	double vxaux,vyaux,vzaux,zaux;
	double Epx,Epy,Epz,Bpx,Bpy,Bpz;
	int n;
	int ppp[ndg][nplim+1],npsinp;
	
	//FFT variables
	double sp1[ndt],sp2[ndt],sp3[ndt],sp4[ndt],raux;
	complex csp1[ndt/2],csp2[ndt/2],csp3[ndt/2],csp4[ndt/2];
	

	
    //Printing parameters in a file
    fprintf(fp0, "General parameters:\nParticle number=%d\nPoints in the spatial grid=%d\nPoints in the temporal grid=%d\n\nCalculation parameters:\nlg=%.3f\tdeltag=%.3f\n\nMaterial parameters:\nxgleft=%.3f\txgright=%.3f\nileft=%d\tiright=%d\n", np, ndg, ndt, lg, deltag, xgleft, xgright, ileft, iright);
	fprintf(fp0,"dt=%e\t Ndt=%d w0=%e\n",deltat,ndt,w0);
	/******************************************************/
	/* Initialize variables
	/* Calculate Densities, Fields and position at the beggining
	/******************************************************/
	

	if (iflagion==0)
	{
		npion=np;
		for(int j=0;j<np;j++)
		{
			qp[j]=-q*n0*(xgright-xgleft)/np;
			mp[j]=mass*n0*(xgright-xgleft)/np;
		}
	}
	else
	{
		npion=0.;
	}
	
	//Writing parameters
	nsc=100;
	nscleft=int(xgleft-0.2*lambda)/deltag;
	nscright=int(xgright+0.2*lambda)/deltag;
	
	for(int j=0;j<nsc;j++)
	{
		idj[j]=nscleft+int((nscright-nscleft)/(nsc-1.)*(j-1.));
	}
	
	
	//parameters of tunnel ionization.
	coef=4.0;
	exponente=2./3.;
	npsinp=0;
	npminimo=10000.;
	cargaminima=n0/100000.;
	
	//Initialize aleatory numbers
	iff=0;
	iseed=1;
	
	//Setting to zero all variables
	vvx=0.;
	vvy=0.;
	vvz=0.;
	
	gvxs=0.;
	gvys=0.;
	gvzs=0.;
	
	aexpt=0.;
	aeypt=0.;
	aezpt=0.;
	
	for(int i=0;i<ndg;i++)
	{
		Ro[i]=0.;
		Rop[i]=0.;
		Jx[i]=0.;
		Jy[i]=0.;
		Jz[i]=0.;
		IntRo[i]=0.;
		IntJx[i]=0.;
		IntJz[i]=0.;
		IPhiP[i]=0.;
		IPhiM[i]=0.;
		
		IAxP[i]=0.;
		IAxM[i]=0.;
		IAyP[i]=0.;
		IAyM[i]=0.;
		Ax[i]=0.;
		Ay[i]=0.;
		Az[i]=0.;
		phi[i]=0.;
		DAx[i]=0.;
		DAy[i]=0.;
		Axold[i]=0.;
		Azold[i]=0.;
		Ex[i]=0.;
		Ey[i]=0.;
		Ez[i]=0.;
		By[i]=0.;
		Bz[i]=0.;
	}
	
	//Initialize particles at rest.	
	if(iflagion==0)
	{
		for(int i=0;i<npion;i++)
		{
			xp[i]=(xgright-xgleft)/npion*1.*i+xgleft;
            vpx[i]=0.;
			vpy[i]=0.;
			vpz[i]=0.;
			double gamma=1./sqrt(1.-(vpx[i]*vpx[i]+vpy[i]*vpy[i]+vpz[i]*vpz[i])/lightC_au/lightC_au);
			vpx[i]=vpx[i]*gamma;
			vpy[i]=vpy[i]*gamma;
			vpz[i]=vpz[i]*gamma;
		}
	}//end if
	
	for(int i=0;i<np;i++)
	{
		deltavx[i]=0.;
		deltavy[i]=0.;
		deltavz[i]=0.;
		deltax[i]=0.;
	}
	
	//Calculate the velocities;
	
 	CalculaRo(xp,np,npion,Ro,ndg,qp, deltag, na, fp3);
	CalculaJ(labelx,xp,vpx,vpy,vpz,np,npion,Jx,Jy,Jz,ndg,qp, deltat, deltag, na, fp2); //labelx is 1 and means 'x'
	CalculaJ(labely,xp,vpx,vpy,vpz,np,npion,Jx,Jy,Jz,ndg,qp, deltat, deltag, na, fp2); //labely is 2 and means 'y'
	CalculaJ(labelz,xp,vpx,vpy,vpz,np,npion,Jx,Jy,Jz,ndg,qp, deltat, deltag, na, fp2); //labelz is 3 and means 'z'
	
	//Density of positive charge
	if(iflagion==0)
	{
		for(int i=0;i<ndg;i++)
			Rop[i]=-Ro[i];
	}
	else
	{
		for(int i=0;i<ndg;i++)
		{
			Rop[i]=0.;
			p0[i]=0.;
			deltap0[i]=0.;
			ppp[i][i]=0.;
		}
		for(int i=ileft;i<iright;i++)
			p0[i]=n0;
		
	}//end if
	
	for(int i=0;i<ndg;i++)
	{
		Axold[i]=Ax[i];
		Ayold[i]=Ay[i];
		Azold[i]=Az[i];
		Exold[i]=Ex[i];
		Eyold[i]=Ey[i];
		Ezold[i]=Ez[i];
		
		Byold[i]=By[i];
		Bzold[i]=Bz[i];
		Roold[i]=Ro[i];
		Ropold[i]=Rop[i];
	}
	
	
	
	/******************************************************/
	/* Central part of the program
	/* Time and particle loop 
	/******************************************************/
	int pacer=0;
	for(int n=0;n<ndt;n++)
	{
		printf("n=%d\n", n);
        pacer+=1;
		t=n*deltat;
		for(int j=1;j<np;j++)
		{
			//inperpola(xp,Ex,Ey,Ez,By,Bz,Epx,Epy,Epz,Bpx,Bpy,Bpz,t,ndg); //Calculate the field in the particle position
			interpola(xp[j],Ex,Ey,Ez,By,Bz,Epx,Epy,Epz,Bpx,Bpy,Bpz,t,ndg,deltag,na);
			double gamma=sqrt(1.+(vpx[j]*vpx[j]+vpy[j]*vpy[j]+vpz[j]*vpz[j])/lightC_au/lightC_au );
			
			aexpt=ae*Epx;
			aeypt=ae*Epy;
			aezpt=ae*Epz;
			
			interpp( &Epy);
			printf("%e ",Epy);
			double sum=Epy;//aexpt+aeypt+aezpt;
			if(sum!=0.)
				printf("+ %e",sum);
			
			//fprintf(fp7,"%e %e %e\n",aexpt,aeypt,aezpt);
			
			abxpt=ae*Bpx;
			abypt=ae*Bpy;
			abzpt=ae*Bpz;
			
			gvxs=vpx[j]+aexpt;
			gvys=vpy[j]+aeypt;
			gvzs=vpz[j]+aezpt;
			
			gamma=sqrt(1.+(gvxs*gvxs+gvys*gvys+gvzs*gvzs)/lightC_au/lightC_au);
		//	fprintf(fp7,"%e\n",gamma);
			dcg=1./lightC_au/gamma;
			
			abxpt=dcg*abxpt;
			abypt=dcg*abypt;
			abzpt=dcg*abzpt;
			
			vvx=gvxs+gvys*abzpt-gvzs*abypt;
			vvy=gvys+gvzs*abxpt-gvxs*abzpt;
			vvz=gvzs+gvxs*abypt-gvys*abxpt;
			
			//Second half rotation and second half acceleration.
			fb=2./(1.+abxpt*abxpt+abypt*abypt+abzpt*abzpt);
			
			abxpt=abxpt*fb;
			abypt=abypt*fb;
			abzpt=abzpt*fb;
			
			gvxs=gvxs+vvy*abzpt-vvz*abypt+aexpt;
			gvys=gvys+vvz*abxpt-vvx*abzpt+aeypt;
			gvzs=gvzs+vvx*abypt-vvy*abxpt+aezpt;
			
			gamma=sqrt(1.+(gvxs*gvxs+gvys*gvys+gvzs*gvzs)/lightC_au/lightC_au);
			
			vpx[j]=gvxs;
			vpy[j]=gvys;
			vpz[j]=gvzs;
			
			xp[j]=xp[j]+deltat*vpx[j]/gamma;
		}//End of the Boris algorithm to solve the equation of motion
		
		//CalculaRo(xp,np,npion,Ro,ndg,qp, deltag, na, fp1);
		CalculaRo_history(xp,np,npion,Ro,ndg,qp, deltag, na, fp1,n);
	    CalculaJ(labelx,xp,vpx,vpy,vpz,np,npion,Jx,Jy,Jz,ndg,qp, deltat, deltag, na, fp2); //labelx is 1 and means 'x'
	    CalculaJ(labely,xp,vpx,vpy,vpz,np,npion,Jx,Jy,Jz,ndg,qp, deltat, deltag, na, fp2); //labely is 2 and means 'y'
        CalculaJ(labelz,xp,vpx,vpy,vpz,np,npion,Jx,Jy,Jz,ndg,qp, deltat, deltag, na, fp2); //labelz is 3 and means 'z'
		
		
		
		/******************************************************/
		/* Calculate Phi, Ax, Ay and Az
		/* By solving the wave equation 
		/******************************************************/


		for(int i=1;i<ndg-1;i++)
		{
			//Time integrals of Ro and J
			
			IntRo[i]=IntRo[i]+(Ro[i]+Rop[i]+Roold[i]+Ropold[i])*0.5*deltat;  //Ros and j are defined in diff time positions
			
			IntJx[i]=IntJx[i]+Jx[i]*deltat;
			IntJy[i]=IntJy[i]+Jy[i]*deltat;
			IntJz[i]=IntJz[i]+Jz[i]*deltat;
			
			phi[i]=2.*pi*lightC_au*(IPhiM[i-1]+IPhiP[i+1]+deltag*IntRo[i] );
			
			Ax[i]=2.*pi*(IAxM[i-1]+IAxP[i+1]+deltag*IntJx[i] );
			Ay[i]=2.*pi*(IAyM[i-1]+IAyP[i+1]+deltag*IntJy[i] );
			Az[i]=2.*pi*(IAzM[i-1]+IAzP[i+1]+deltag*IntJz[i] );
		}
		
		IntRo[0]=IntRo[0]+(Ro[0]+Rop[0]+Roold[0]+Ropold[0])*0.5*deltat;
		
		IntJx[0]=IntJx[0]+Jx[0]*deltat;
		IntJy[0]=IntJy[0]+Jx[0]*deltat;
		IntJz[0]=IntJz[0]+Jx[0]*deltat;
		
		phi[0]=2.*pi*lightC_au*(IPhiP[1]+deltag*IntRo[0]);
		Ax[0]=2.*pi*(IAxP[1]+deltag*IntJx[0]);
		Ay[0]=2.*pi*(IAyP[1]+deltag*IntJy[0]);
		Az[0]=2.*pi*(IAzP[1]+deltag*IntJz[0]);
		
		IntJx[ndg-1]=IntJx[ndg-1]+Jx[ndg-1]*deltat;
		IntJy[ndg-1]=IntJy[ndg-1]+Jx[ndg-1]*deltat;
		IntJz[ndg-1]=IntJz[ndg-1]+Jx[ndg-1]*deltat;
		
		phi[ndg-1]=2.*pi*lightC_au*(IPhiP[ndg-2]+deltag*IntRo[ndg-1]);
		Ax[ndg-1]=2.*pi*(IAxP[ndg-2]+deltag*IntJx[ndg-1]);
		Ay[ndg-1]=2.*pi*(IAyP[ndg-2]+deltag*IntJy[ndg-1]);
		Az[ndg-1]=2.*pi*(IAzP[ndg-2]+deltag*IntJz[ndg-1]);
		
		for (int i=0;i<ndg;i++)
		{
			DAx[i]=(Ax[i]-Axold[i]);
			DAy[i]=(Ay[i]-Ayold[i]);
			DAz[i]=(Az[i]-Azold[i]);			
		}
		
		//Spatial integrals.
		for (int i=0;i<ndg-1;i++)
		{
			IPhiP[i]=IPhiP[i+1]+IntRo[i]*deltag;
			IAxP[i]=IAxP[i+1]+IntJx[i]*deltag;
			IAyP[i]=IAyP[i+1]+IntJy[i]*deltag;
			IAzP[i]=IAzP[i+1]+IntJz[i]*deltag;
		}
		
		IPhiP[ndg-1]=IntRo[ndg-1]*deltag;
		IAxP[ndg-1]=IntJx[ndg-1]*deltag;
		IAyP[ndg-1]=IntJy[ndg-1]*deltag;
		IAzP[ndg-1]=IntJz[ndg-1]*deltag;
		
		for(int i=ndg-1;i>=1;i--)
		{
			IPhiM[i]=IPhiM[i-1]+IntRo[i]*deltag;
			IAxM[i]=IAxM[i-1]+IntJx[i]*deltat;
			IAyM[i]=IAyM[i-1]+IntJy[i]*deltat;
			IAzM[i]=IAzM[i-1]+IntJz[i]*deltat;
		}

		IPhiM[0]=IntRo[0]*deltag;
		IAxM[0]=IntJx[0]*deltat;
		IAyM[0]=IntJy[0]*deltat;
		IAzM[0]=IntJz[0]*deltat;
		
		for(int i=0;i<ndg;i++)
		{
			Axold[i]=Ax[i];
			Ayold[i]=Ay[i];
			Azold[i]=Az[i];

			Roold[i]=Ro[i];
			Ropold[i]=Rop[i];
		}
		
		/******************************************************/
		/* End of the
		 /* Calculation Phi, Ax, Ay and Az
		 /* By solving the wave equation 
		 /******************************************************/
		
		
//		CalculaEx(Ex,phi,DAx,ndg,t+deltat);		
		CalculaEx(Ex,phi,DAx,ndg,t+deltat,deltag,deltat);
//		CalculaEy(Ex,DAy,ndg,t+deltat);
		CalculaEy(Ey,DAy,ndg,t,deltag,deltat,xg,xgleft,wp,kp,E0,k0,w0);
//		CalculaEz(Ez,DAz,ndg,t+deltat);
		CalculaEz(Ez,DAz,ndg,t+deltat,deltag,deltat);		
//		CalculaBy(By,Az,ndg,t+deltat);
		CalculaBy(By,Az,ndg,t+deltat,deltag,deltat);
//		CalculaBz(Bz,Ay,ndg,t+deltat);
		CalculaBz(Bz,Ay,ndg,t+deltat,deltag,deltat,xg,xgleft,wp,kp,E0,k0,w0);			
		
		
		
		
		
		if(iflagion!=0)
		{
			//Ionization();
			//Calculap0();
		}
		
		
		//Calculate the transmitted and reflected field.
		sp1[n+1]=Ey[int((xgleft-lambda)/deltag)]-E0y(xgleft-lambda,t+deltat,xgleft,wp,kp,E0,k0,w0);
		sp2[n+1]=Ey[int((xgright+lambda)/deltag)];
		
		for(int i=1;i<=ndg/10;i++)
		{
			fprintf(fp6,"%e\n",Ey[i*10]);
		}
		
	}//end time loop
	
	//Printing the transmited field
	for(int i=1;i<=ndt;i++)
    {
		fprintf(fp6,"%d\t%e\n",i, sp2[i]);

    }
	
	
	
	fclose(fp0);
	fclose(fp1);
    fclose(fp2);	
	fclose(fp3);
    fclose(fp4);	
	fclose(fp5);
    fclose(fp6);		
	
	
}


void interpp(double *a)
{
	*a=100.;
}
